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Recent experimental observations have demonstrated interesting instability phenomenon during 
thermal drawing of microstructured glass/polymer fibers, and these observations motivate us to 
examine surface-tension-driven instabilities in concentric cylindrical shells of viscous fluids. In this 
paper, we focus on a single instability mechanism: classical capillary instabilities in the form of radial 
fluctuations, solving the full Navier-Stokes equations numerically. In equal- viscosity cases where an 
analytical linear theory is available, we compare to the full numerical solution and delineate the 
regime in which the linear theory is valid. We also consider unequal-viscosity situations (similar 
to experiments) in which there is no published linear theory, and explain the numerical results 
with a simple asymptotic analysis. These results are then applied to experimental thermal drawing 
systems. We show that the observed instabilities are consistent with radial-fluctuation analysis, but 
cannot be predicted by radial fluctuations alone — an additional mechanism is required. We show 
how radial fluctuations alone, however, can be used to analyze various candidate material systems 
for thermal drawing, clearly ruling out some possibilities while suggesting others that have not yet 
been considered in experiments. 



I. INTRODUCTION 

The classical capillary instability, the breakup of a 
cylindrical liquid thread into a series of droplets, is per- 
haps one of the most ubiquitous fluid instabilities and 
appears in a host of daily phenomena [Ij ^2, from glass- 
wine tearing and faucet dripping to ink-jet print- 
ing. The study of capillary instability has a long his- 
tory. In 1849, Plateau attributed the mechanism to sur- 
face tension: the breakup process reduces the surface en- 
ergy [4 . Lord Rayleigh pioneered the application of lin- 
ear stability analysis to quantitatively characterize the 
growth rate at the onset of instability, and found that a 
small disturbance is magnified exponentially with time 
[S]. Subsequently, Tomotika investigated the effect of 
viscosity of the surrounding fluid, showing that it acts 
as a deterrent to slow down the instability growth rate 
[B]. Stone and Brenner investigated instabilities in a 
two-fluid cylindrical shell geometry with equal viscosities 
[7 . Many additional phenomena have been investigated, 
such as the cascade structure in a drop falling from a 
faucet [8J , steady capillary jets of sub-micrometer diam- 
eters |9j, and double-cone neck shapes in nanojets |10| . 
Capillary instability also offers a means of controlling and 
synthesizing diverse morphological configurations. Ex- 
amples include: a long chain of nanospheres generated 
from tens-of-nanometer diameter wires at the melt state 
[111 I12| ; polymer nanorods formed by annealing poly- 
mer nanotubes above the glass transition point |T3]; and 
nanoparticle chains encapsulated in nanotubes generated 
by reduction of nanowires at a sufficiently high temper- 
ature [141. Moreover, instabilities of fluid jets have nu- 
merous chemical and biological applications f l5l[TB] . (An 
entirely different instability mechanism has attracted re- 
cent interest in elastic or visco-elastic media, in which 
thin sheets under tension form wrinkles driven by elastic 
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Recently, a new class of glass and polymer microstruc- 
tured fibers, which are characterized by an embedded 
geometry of concentric cylindrical shells, has emerged 
for a variety of applications in optics |20fl22| . Uniform- 
thickness cylindrical shells (made of glass materials such 
as As2Se3 and AS2S3), down to sub- micrometer or even 
nanometer thicknesses, have been successfully fabricated 
in glass materials, by a drawing process whereby a large- 
scale "preform" is heated and pulled into a long thread, 
as depicted in Fig. [T] On the other hand, as the shell 
thickness is further reduced towards the nanoscale, the 
thin cylindrical shell (made of the glass material sele- 
nium, Se), is observed to break up into an ordered array 
of filaments; that is, the breakup of the cylindrical shell 
occurs along the azimuthal direction while the axial con- 
tinuity remains intact [2TJ |22| . 

Fiber drawing of cylindrical shells and other mi- 
crostructured geometries clearly opens up rich new ar- 
eas for fluid instabilities and other phenomena, and it is 
necessary to understand these phenomena in order to de- 
termine what structures are attainable in drawn flbers. 
Although it is a striking example, the observed azimuthal 
breakup process appears to be a complicated matter — 
because surface tension does not produce azimuthal in- 
stability in cylinders |23l |21] or cylindrical shells |25] . 
the azimuthal breakup must be driven by the fiber draw- 
down process and / or by additional physics such as visco- 
elastic effects or thermal gradients. Simulating the en- 
tire draw-down process directly is very challenging, how- 
ever, because the lengthscales vary by 5 orders of mag- 
nitude from the preform (cm) to the drawn layers (sub- 
/im). Another potentially important breakup process, 
one that is amenable to study even in the simplified case 
of a cylindrical-shell geometry, is axial instability. Not 
only does the possibility of axial breakup impose some 
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limits on the practical materials for fiber drawing, but 
understanding such breakup is arguably a prerequisite 
to understanding the azimuthal breakup process, for two 
reasons. First, the azimuthal breakup process produces 
cylindrical filaments, and it is important to understand 
why these filaments do not exhibit further breakup into 
droplets (or under what circumstances this should oc- 
cur). Second, it is possible that the draw-down process 
or other effects might couple fiuctuations in the axial and 
azimuthal directions, so understanding the timescales of 
the axial breakup process is necessary as a first step in 
evaluating whether it plays any physical role in driving 
other instabilities. 

Therefore, as a first step towards understanding the 
various instability phenomena in drawn microstructured 
(cylindrical shell) fibers, we investigate the impact of a 
single mechanism: classical Plateau-Rayleigh instability 
leading to radial fiuctuations in cylindrical fiuids, here in- 
cluding the previously unstudied case of concentric shells 
with different viscosity. To isolate this mechanism from 
other forms of instability, we consider only cylindrically 
symmetrical geometries, which also greatly simplifies the 
problem into two dimensions (r, z). (Indeed, this model 
is a satisfactory description of the cylindrical filaments.) 
We model this situation by direct numerical simulation 
via a finite-element method for the Navier-Stokes (NS) 
equations. After first reproducing the known results from 
the linear theory for equal viscosity, we are able to explore 
the regime beyond the limits of the linear theory, which 
we show is accurate up to fiuctuations of about 10% in ra- 
dius. For the case of unequal- viscosity shells, where there 
is not as yet a linear-theory analysis, we show that the 
instability timescale interpolates smoothly between two 
limits that can be understood by dimensional analysis. 
Finally, we apply our results to the experimental fiber- 
drawing situation, where we obtain a necessary (but not 
sufficient) condition for stability that can be used to guide 
the materials selection and the design of the fabrication 
process by excluding certain materials combinations from 
consideration. Our stability criterion is shown to be con- 
sistent with the experimental observations. We also find 
that the stability of the resulting filaments is consistent 
with the Rayleigh-Tomotika model. 

Motivated by the desire to understand the range of at- 
tainable structures and improve their performance, pre- 
vious theoretical study of microstructured fiber drawing 
has considered a variety of situations different from the 
one considered here. Air-hole deformation and collapse 
was explored by numerical analysis of the continuous 
drawing process of microstructured optical fibers [26 , 27J. 
Surface-tension effects have been studied for their role 
in determining the surface smoothness and the result- 
ing optical loss of the final fibers The modeling of 
non-circular fibers has also been investigated in order to 
design fiber-draw process for unusual fiber shapes with a 
square or rectangular cross-section (29]. 

This paper is organized as follows. We provide more 
background on microstructured fibers in section |IT] and 



review the governing equations and dimensionless groups 
pertaining to fiber thermal-drawing processing in sec- 
tion lllll Section HVl describes the numerical finite-element 
approach to solving the Navier-Stokes equations, and 
section V presents our simulation results for cylindrical 
shell. Section VI presents a radial stability map to guide 
materials selection, which is established by linear-theory 
calculations dependent on the shell radius, thickness, and 
viscosities. In section |VII[ we discuss the applications 
of capillary instability to our microstructured fibers for 
viscous materials selection and limits on the ultimate 
feature sizes. In particular, section |VIID| gives a sim- 
ple geometrical argument to explain why any instability 
mechanism in fiber drawing will tend to favor azimuthal 
breakup (into filaments) over axial breakup (into rings or 
droplets) . 



II. FEATURE SIZE IN COMPOSITE 
MICROSTRUCTURED FIBERS 

Most optical fibers are mainly made of a single mate- 
rial, silica glass. Recent work, however, has generalized 
optical fiber manufacturing to include microstructured 
fibers that combine multiple distinct materials includ- 
ing metals, semiconductors, and insulators, which expand 
fiber-device functionalities while retaining the simplicity 
of the thermal-drawing fabrication approach [20 . For 
example, a periodic cylindrical-shell multilayer structure 
has been incorporated into a fiber to guide light in a hol- 
low core with significantly reduced loss for laser surgery 
|30| . The basic element of such a fiber is a cylindrical 
shell of one material (a chalcogenide glass) surrounded by 
a "cladding" of another material (a thermoplastic poly- 
mer), as depicted in Fig. [T| The fabrication process has 
four main steps to create this geometry, (i) A glass film is 
thermally evaporated onto a polymer substrate, (ii) The 
glass/polymer bi-layer film is tightly wrapped around a 
polymer core, (iii) Additional layers of protective poly- 
mer cladding are then rolled around the structure, (iv) 
The resulting centimeter-diameter preform is fused into a 
single solid structure by heating under vacuum. The solid 
preform is then heated into a viscous state and stretched 
into an extended millimeter-diameter fiber by the appli- 
cation of axial tension, as shown in Figjl] 

Uniform layer thicknesses, down to micrometers or 
even to nanometers, have been successfully achieved in 
fibers by this method. Mechanically fiexible fibers with 
a uniform diameter have been produced, as shown in 
Fig. ^a). Fig. |2jb) shows a typical Scanning Electron 
Microscope (SEM) micrograph of a cross-section with 1 
mm fiber diameter. Magnified SEM micrographs in Fig. 
|2jc) and (d), for two different fibers, reveal high-quality 
multiple-layer structures with thicknesses on the order of 
micrometers and tens of nanometers, respectively. 

Ideally, the cross-section of the resulting fibers retains 
the same structure and relative sizes of the components as 
in the preform. If the layer thickness becomes too small. 
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(thermal heating)' 




Figure 1: Optical-fiber thermal drawing. Preform is heated 
at elevated temperature to viscous fluid, and stretched into 
extended flbers by applied tension. This preform is specially 
designed with a thin cylindrical shell in polymer matrix. 




Figure 2: SEM micrographs of cylindrical shells in fiber, (a) 
Photograph of fiber, (b) SEM of fiber cross-section. Magni- 
fied view of multilayer structures reveals the thickness of mi- 
crometer (c) and tens of nanometers (d) , respectively. Bright 
and dark color for glass and polymer in SEM, respectively, (e) 
showing layer breakup into circles in the fiber cross-section, 
while (f) presenting the continuous filaments obtained from 
fiber after dissolving polymer matrix. 



however, we observed the layers to break up azimuthally, 
as shown in the cross-section of Fig. [2f^e), while conti- 
nuity along the axial direction remains intact. In this 
fashion, a thin cylindrical shell breaks into filament ar- 
rays embedded in a fiber [2T1 [55]. After dissolving the 
polymer cladding, the resulting separated glass filaments 
are shown in Fig. [2jf) [21 . 

As a first step towards understanding these observa- 
tions, we will consider simplified geometries consisting of 
a thin cylindrical shell in the cladding matrix. The ques- 



tion, here, is whether the classical capillary instability 
from radial fluctuations is relevant to the attained thin 
uniform shells and the observed azimuthal breakup at 
further reduced thicknesses. Moreover, we want to know 
whether classical capillary instability can provide guid- 
ance for materials selections and fabrication processing 
in our microstructured fibers. 



III. GOVERNING EQUATIONS 

During thermal drawing, the temperature is set above 
the softening points of all the materials, which conse- 
quently are in the viscous fluid state to enable polymer 
and glass codrawing. To describe this fluid flow, we con- 
sider the incompressible Navier-Stokes (NS) equations 
ED: 



p [dtu + {u • V) w] = 



■q- 



(1) 

where u is velocity, p is pressure, p is materials density, 
T] is viscosity, 7 is interfacial tension between glass and 
polymer, S is the delta function (S — 1 at the polymer 
and glass interface and S = otherwise), k is curvature 
of interface, and n is a unit vector at interface. 

To identify the operating regime of the fiber drawing 
process, we consider the relevant dimensionless numbers. 
The Reynolds number (Re), Froude number (Fr), and 
capillary number(Ca) are: 



Re = 



pUh 



— ,Fr = — ,Ca = — , 
77 gh 7 



(2) 



where p w 10"^ kg/m'^ is the density of the materials, g w 
10 m/s^ is gravity, U ~ 5 mm/s is drawing speed, 77 « 10^ 
Pa • s is viscosity, /i « 100 nm is the layer thickness, and 
7 = 0.1 N/m is surface tension between polymer and 
glass |32| 133] . Therefore, these dimensionless numbers 
in a typical fiber draw are Re « 10"^*^, Fr « 10^ and 
Ca w 10^. Small Re number, large Fr number, and large 
Ca number imply a weak inertia term, negligible gravity, 
and dominant viscosity effects, respectively. In addition, 
since the fiber diameter is Z) « 1mm and the length of the 
neck-down region is L « 10 cm, the ratio D/L sa 1/100 
is much less than 1, and thus the complicated profile of 
neck-down cone is simplified into a cylindrical shape for 
the purpose of easier analysis. 



IV. SIMULATION ALGORITHM 

To develop a quantitative understanding of capillary 
instability in a cylindrical-shell geometry, direct numer- 
ical simulation is performed using the finite element 
method. In order to isolate the effect of radial fluctua- 
tions, we impose cylindrical symmetry, so the numerical 
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simulation simplifies into a 2D problem in the (r, z) plane 
(where the V operations are replaced by their cylindri- 
cal forms). Although the low Re number means that 
one could accurately neglect the (u ■ V)u nonlinear iner- 
tia term and solve only the time-dependent Stokes prob- 
lem, we solve the full NS equations because of the conve- 
nience of using available software supporting these equa- 
tions and also to retain the flexibility to consider large 
Re regimes in the future. 

Our simulation algorithm is briefly presented in Fig. 
[Sj^a). The geometry of a cylindrical shell is defined by 
two interfaces I and II in Fig. |3jb). The flow fleld {u,p) 
is calculated from Navier-Stokes equations, as seen in 
Fig. |3jc). Consequently, this new flow field generates 
interface motion, resulting in an updated interface and 
an updated flow field. By these numerical iterations, the 
interface gradually evolves with time. 

More specifically, the time-dependent interfaces I and 
II can be expressed as follows: 

ri{z,t) = ro + ei{z,t), 

where ro, 2ro are the unperturbed radii of interface I and 
II, ei(z, t), £2(2;, i) denote the interfacial perturbations 
that grow with time. 

The initial perturbations at interface I and II are given 
by cosine- wave shapes: 

ei(z,0)= eoCOs(27rz/A), 
£2(2,0) = 2eocos(27rz/A), 

where eo,2eo are the perturbation amplitudes, and A is 
the initial perturbation wavelength. 

Numerical challenges in the simulations arise from the 
nonlinearity, moving interfaces, interface singularities, 
and the complex curvature |TJ [33]. A level-set func- 
tion 4){x^ t) is coupled with the NS equations to track the 
interface (see Appendix 1) |35ff57] . where the interface 
is located at the ^ = 0.5 contour and the evolution is 
given by u via: 



(a) Interface motion 



' 1 Navier - Stol(es equation i 

Level Set J Flow field 




Figure 3: Simulation algorithm, (a) Schematic of algorithm, 
(b) Interfaces I and II of cylindrical shell are defined by a level 
set function, (c) Flow field is determined by the NS equations 
(color scale for pressure and arrows for fluid velocity). 



V. SIMULATION RESULTS 

In the low Re regime, the NS equations are linear with 
respect to the fiow velocity u because the inertia term 
{u-V)u is negligible, but they are not linear with respect 
to the geometry of the interface shape except in the limit 
of small perturbations. A linear theory for small geo- 
metric perturbations was formulated for the cylinder by 
Tomotika fG] and for a cylindrical shell with a cladding 
of equal viscosity by Stone and Brenner [7]. Our full 
simulation should, of course, reproduce the results of the 
linear theory for small perturbations, but also allows us 
to investigate larger perturbations beyond the regime of 
the linear theory, and we can also consider the more gen- 
eral case of a cylindrical shell with a cladding of unequal 
viscosity. 



A. Numerical convergence 



0t -(- u • V(/) = 0. 



(5) 



The local curvature (k) at an interface is given in terms 
of (f) by: 



V0 



|V0| 



_ {4>l(t)zz + 4'l<Prr) + [4>l + 4'l)<t)r/r - <PAz{<Prz + 4>zr) 
(02+^2)3/2 

(6) 

For numerical stability reasons, we add an artificial diffu- 
sion term proportional to a small parameter D in Eq. [s] 
(see Appendix 4) ; the convergence as 13 is discussed 
in the next section. 



In order to establish the accuracy of the simulation, 
we first investigated the convergence with respect to the 
various approximations: the time integration tolerance 
(— >■ 0), the mesh resolution (— >■ 00), and an artificial dif- 
fusion term D (-> 0). The time integration is performed 
by a fifth-order backward differentiation formula, charac- 
terized by specified relative and absolute error tolerances 
(errrci = lO^** and errabs = 10^^) with acceptable accu- 
racy (w 10^'^). Fig. [4 a) shows numerical convergence 
of the instabiltiy amplitude mh\z £i{z^t) at time t=300 
seocnds as a fucntion of the number of triangular mesh 
elements. From last four data points in (a), the con- 
vergence of |Ae|/ro with respect to the mesh resolution, 
compared with a very fine mesh of 20-10'* mesh elements, 
is presented in Fig. |4f|b). Because of the PDE's discon- 
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Figure 4: Numerical convergence with respect to various 
approximations in the simulation, (a) Numerical conver- 
gence shown by instability amplitude at time t=300 seoc- 
nds vs number of mesh elements; (b) established errors (|Ae| 
corresponding to higher resolution) versus resolution (solid 
line=expected asymptotic rate); (c) the artificial diffusion 
term for numerical stability, and (d) the numerical errors rep- 
resented by deviations from volume conservation. 




20 40 60 80 100 

Time (units of rpiyy) 



tinuous coefficient, we expect Ist-order convergence: that 
is error '-^ 1/ (mesh resolution) , corresponding to slope 
— 1 (solid line) in Fig. |4|b): the data points appear to 
be approaching this slope, but we are not able to mea- 
sure errors at high enough resolution to verify the con- 
vergence rate conclusively. However, we found 10 • 10** 
mesh elements, corresponding to a typical element diam- 
eter « O.Olro , yield good accuracy (« 10^"^ error). The 
physical results are recovered in the limit Z? — 0, and 
so we must establish the convergence of the results as D 
is reduced and identify a D small enough to yield accu- 
rate results while still maintaining stability. Fig. |4|^c) 
demonstrates this convergence in the simulation by cal- 
culating the amplitude difference jAej/ro compared with 
D = 0.5i:»o (-Do = lO^i'^mVs in Appendix 1), and 3 dig- 
its of accuracy are attained by I? = Dq. Accordingly, the 
numerical parameters throughout the following simula- 
tions are errrei = 10""*, errahs = 10~^ , D = I.ODq, and 
10 • lO"' mesh elements. As an additional check, the vol- 
ume of the incompressible fluid shell between the two in- 
terfaces should be conserved. As shown in Fig. |4|^d), the 
numerical errors in the simulation leads to only a small 
volume deviation, {\V{t) - ViO)\/V{0) = \AV{t)\/V{0)), 
of less than 0.4%. We conclude that the finite-element 
level-set approach is numerically reliable for the simula- 
tion of capillary instabilities. 



Figure 5: (a) Snapshot of the flow field and interfaces during 
instability evolution. Color scale for pressure, arrows for fiuid 
velocity, (b) The instability amplitude grows exponentially 
with time at interfaces I and II. 



B. Instability evolution 

The evolution of a capillary instability can be ob- 
tained from direct numerical simulation. Fig. [sj^a) (i)-(v) 
presents snapshots of the flow field and interface, (i) Ini- 
tially, the pressure of the inner fluid is higher than that 
of the outer fluid due to Laplace pressure (p = ^k) origi- 
nating from azimuthal curvature of the cylindrical geom- 
etry at interfaces I and II. (ii)-(iv) The interfacial per- 
turbations generate an axial pressure gradient Ap, and 
hence a fluid flow occurs that moves from a smaller-radius 
to a larger-radius region for the inner fluid. Gradually 
the amplitude of the perturbation is amplified, (v) The 
shrunk smaller-radius and expanded larger-radius regions 
of inner fluid further enhance the axial pressure gradient 
Ap, resulting in a larger amplitude of the perturbation. 
As a result, the small perturbation is exponentially am- 
plifled by the axial pressure gradient. 

Fig. |5jb) shows the instability amplitude miriz e{z,t) 
at interfaces I and II growing with time exponentially on 
a semilog scale in the plot. The instability time scale in 
the simulation is about 372 ± 3 sec by fltting the curve 
to an exponential e ~ e*/"^. Linear theory predicts time 
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Figure 6: Growth of perturbation amplitude as a function 
of time. At small amplitudes below 10% r, the perturba- 
tion grows exponentially as expected from linear theory. At 
large amplitude above 10% r, significant deviations from lin- 
ear theory occur. Scale invariance, with various values of r 
superimposed, is expected in the Stokes regime with low Re 
number. 

scale of r w 334 sec |7] (see Appendix 2). Thus, simula- 
tion is consistent with the linear theory, in terms of the 
instability time scale. 



C. Beyond the linear theory 

More quantitatively, for small perturbations or short 
times, exponential growth of the perturbation amplitude 
is expected from the linear theory. However, at later 
times, when the perturbation has grown sufficiently, one 
expects the linear theory to break down. We observe sig- 
nificant deviations from the linear theory if the amplitude 
is above 10% of the cylinder radius r (corresponding to 
time 40rr//7), as shown in Fig. [6] (taking interface I as 
an example). For example, at perturbation amplitudes 
of around 25% of r (time lOrij/j), the deviations of the 
simulation data from the linear theory is almost 20% . 

The length scaling of the instability is studied as well. 
After rescaling time and distance by the same factor (so 
that the initial radius is scaled by 0.25,0.5, 1,2, and 4), 
all the data collapse onto a single master curve, as shown 
in Fig. [6] This length scaling originates from the low Re 
regime, which implies that the NS equations are approxi- 
mately linear with respect to u; this gives rise to the well 
known scale invariance in the Stokes regime [31j . 



D. Unequal viscosities 

We now investigate the dependence of instability 
timescale on cladding viscosity (r/ciad) with a fixed shell 



Figure 7: Time-dependent perturbation amplitude curves for 
various cladding viscosities (77ciad), but with a fixed shell vis- 
cosity (r/shcii = lO^Pa- s). 



viscosity (77shcii), in order to help us to identify suit- 
able cladding materials for fiber fabrication. As viscosity 
slows down fluid motion, the low- viscosity cladding has 
a faster instability growth rate and a shorter time scale, 
while the high-viscosity cladding has a slower instability 
growth rate and a longer instability time scale. Fig. [7] 
shows the time-dependent perturbation amplitude curves 
for various viscosity contrast 77ciad/?7shcii by changing the 
cladding viscosity (?7shoU — 10^ Pa-s). Instability time 
scale for the each given viscosity contrast is obtained by 
exponentially fitting the curves in Fig. [7j Instability time 
scale (r) as a function of viscosity contrast is presented 
in Fig. Ill 

The existing linear theory has only been solved in case 
of equal viscosity, and predicts that the instability time 
scale is proportional to the viscosity r ^ 77 [7, . We obtain 
a more general picture of the instability time scale for un- 
equal viscosity by considering two limits. In the limit of 
negligible cladding viscosity, ?7ciad — ^ 0, the instability 
time scale should be determined by rysheii, and from di- 
mensional analysis should be proportional to rrjshcii/"/ 1 
assuming that the inner and outer radius are compara- 
ble and so we take r to be the average radius. In the 
opposite limit of ?7ciad ~^ 00 , the time scale should be 
determined by rydad and hence should be proportional to 
fVciad/j- In between these two limits, we expect the time 
scale to smoothly interpolate between the rryshoii/7 and 
f"nciad/l scales. Precisely this behavior is observed in 
our numerical calculations, as shown in Fig. [8j 



VI. ESTIMATE OF RADIAL INSTABILITY 
TIMESCALE 

In the previous section, we surveyed the capillary in- 
stability of a concentric-cylindrical shell by numerical 
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Figure 8: An instability time scale (r) for unequal viscosity 
with a fixed shell viscosity (77shGii)- In the limit of 77ciad — )■ 0, r 
is determined by ryshcii and approaches to a constant. In the 
opposite limit of rydad — >■ oo, r should be determined by T^dad 
and is linearly proportional to T^ciad • Between these two limits 
of ?7ciad, T smoothly interpolates between the corresponding 
time scales. 



simulations. We proceed to apply these results to in- 
vestigate the instability time scale dependent on various 
physical parameters including geometry (e.g., radius and 
shell thickness) and materials properties (viscosity) in a 
broad range, exploiting the accuracy of simple dimen- 
sional analysis demonstrated in the previous section and 
providing guidance of our fiber drawing. 

The calculated instability time scale (r) for different 
values of the radius (r) and the viscosity {rf) is displayed 
in Fig. [9] The cross-sectional geometry in the calcu- 
lation is shown in the inset: interface I is located at 
radius r, the cylindrical-shell thickness is h, and inter- 
face II is at radius R = r + h. The interfacial tension in 
the calculations was set to 7 = O.lN/m, which was the 
measured interfacial tension between thermoplastic poly- 
mer and chalcogenide glass used in our microstructured 
fibers |32]. 

Two cases are considered: one is equal viscosity 
(??sheii = ''?ciad), and the other is unequal viscosity 
(»?sheii 7^ Vcind)- In the case of Tyshdi = ?7ciad, the in- 
stability time scale is calculated exactly from Stone and 
Brenner's linear theory 



rjr 



7maxA*(A,i?/r)' 



(7) 



where the fastest growth factor max^ ^'(A, i?/r) was 
found by searching numerically within a wide range of 
wavelengths A for a certain value oi R/r (Eq. 23 in Ap- 
pendix 2). Fig. [9] plots this time scale versus radius for 
T] = lO'^Pa • s corresponding to As2Se3-PES, compared to 
the dwelling time Towelling ~ 100 sec which is defined by 
the time of materials in viscous state before exiting hot 



Figure 9: Radial stability map. Linear theory calculations 
of the instability time scale (r), which is dependent on the 
radius, thickness, and viscosity. Inset shows cross-sectional 
geometry of cylindrical shell. In our experiments, the dwelling 
time of thermal drawing is around rdwciiing ~ 100 sec, and 
the fiber radius r « 500 /im. Radially stable region is shaded 
yellow for r > Tdwoiiing, while the unstable region corresponds 

to r < Tdwclling. 



furnace to be frozen in fiber during thermal drawing |33] . 

In the other case of r]shcU 7^ ?/ciad (in the 7?ciad/?ysheU > 
1 regime as discussed in Section V D I , the instability time 



scale can be roughly estimated from dimensional analy- 
sis. Although dimensionless analysis does not give the 
constant factor, for specificity, we choose the constant 
coefficient from the Tomotika model ||6j. 



7maxA [(1 - x^)^{x,riciad/'nsheu)] ' 



(8) 



where the fastest growth factor maxA[(l — 
x'^)^{x,r]ciadhsheu)] was found numerically by 
searching a wide range of wavelengths {x — 2T:r/X) 
[$(x, Tyciad/'yshcu) IS a Complicated implicit function of 
wavelength and viscosity contrast given in Ref. pi]. Fig. 
[9] plots the time scale for Tysheii = 10 Pa • s, ?7ciad = 10^ 
Pa • s, corresponding to Se-PSU showing that the ob- 
served stability of shells of radius w 250/xm is consistent 
with the radial stability criterion (r > Tdwciiing)- On the 
other hand, if ryciad is reduced to 10^ Pa • s with the same 
shell materials, corresponding to Se-PE, we predict that 
radial ffuctuations alone will render the shell unstable 
for any radius r < 1 mm. 



VII. APPLICATIONS IN MICROSTRUCTURED 
FIBERS 

In this section, we consider in more detail the appli- 
cation of these analyses to understand observed experi- 
mental results, and in particular the observed stability 
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(or instability) of thin shells and filaments. We exam- 
ine whether our stability analysis can provide guidance 
in materials selection and in the the understanding of 
attainable feature sizes. Because we only considered ra- 
dial fluctuations, our analysis provides a necessary but 
not sufficient criterion for stability. Therefore, the rel- 
evant questions are whether the criterion is consistent 
with observed stable structures, whether it is sufficient to 
explain the observed azimuthal breakup, and what mate- 
rials combinations are excluded. Below, in Sec | VII A] we 
consider the application of radial stability analysis to the 
observed stability or instability of cylindrical shells. In 
Sec I VII 5] we look at the impact on materials selections. 
Finally, in Sec | VII C| we show that the observed stability 
of the resulting nanoscale filaments is consistent with the 
Tomotika model. 



Potentially suitable 



A 

ASjSj-PEI 




10 



10° 10' 10' 

Shell viscosity (Pa s) 



10' 



A. Comparison with observations for cylindrical 
shells 



The radial stability map of Fig. |9]is consistent with the 
experimental observations (the default cylindrical shell 
radius « 250^m). First, the map predicts that feature 
sizes down to submicrometers and hundreds of nanome- 
ters are consistent with radial stability for the equal- 
viscosity materials combination of 77ciad — VshcU ~ 10^ 
Pa • s, which corresponds to As2Se3-PES or AS2S3-PEI. 
For As2Se3-PES, Figj2] (c) shows that a shell thickness 
of AS2S3 of 1 fim is obtained; in other work, layers of 
AS2S3 down to 15 nm have been achieved as well |21| . 
For AS2S3-PEI, Figj2] (d) demonstrates a thickness of 
AS2S3 down to 32 nm. Second, the map is consistent with 
thicknesses down to submicrometers for unequal- viscosity 
materials with rydad — 10^ Pa • s, jysheii — 10 Pa • s, which 
corresponds to Se-PSU. Se layers with thickness on the 
order of 1 /im have been demonstrated in Se-PSU fiber 

m- 

The radial stability map, nevertheless, is not sufficient 
to explain the azimuthal instability at further reduced 
thicknesses down to tens of nanometers. The stability 
map of Fig. [9] predicts that a Se layer in a Se-PSU 
combination should be radially stable down to tens of 
nanometers. However, we found in experiments [Fig. [2]Je, 
f)] that a Se shell with thickness < 100 nm breaks up into 
continuous filament arrays |21l [22] , which means that the 
mechanism for this breakup is distinct from that of purely 
radial fiuctuations. For another As2Se3-PES materials 
combination, this filamentation of As2Se3 film was also 
observed as the thickness is reduced down to 10 nm |21| . 
Future work will elucidate this filamentation mechanism 
by performing 3D numerical simulation to explore the 
azimuthal fiuctuations. 



Figure 10: Calculated shell-cladding viscous materials selec- 
tion map during thermal drawing (rdwciiing = 100 sec). A red 
line for instability time for dwelling time r = rdwciung. The 
shaded region above the red line indicates potentially suit- 
able viscous materials combination (r > rdwoiiing), those in 
which radial instabilities alone do not cause breakup. The 
region below the red line indicates radially unstable materi- 
als combinations (r < Tdwciiing), such as Se — PE materials 
combination, which are unstable for thermal drawing. 



B. Materials selection 

Given the viscosities and surface tension of a particu- 
lar material pair, we can use Fig. |9] to help determine 
whether that pair is suitable for drawing: if it is radially 
unstable, then it is almost certainly unsuitable (unless 
the process is altered to somehow compensate) , whereas 
if it is radially stable then the pair is at least poten- 
tially suitable (if there are no other instabilities). Vis- 
cosities of materials used in fiber drawing are obtained 
as follows: the viscosities of semiconductor glasses (Se, 
As2Se3, AS2S3) are calculated from an empirical Arrhe- 
nius formula at the associated temperature during ther- 
mal drawing (more details in Appendix 3); several ther- 
moplastic polymers (PSU, PES, and PEI) have similar 
viscosities Typoiymer ~ 10^ Pa • s during fiber drawing |33| ; 
and the viscosity of the polymer PE is 10'^ Pa • s at tem- 
perature T = 250° C [33. The polymer-glass surface ten- 
sion is typically 7 = 0.1 N/m [32] for all of these materi- 
als. Assuming a cylindrical shell of radius « 250 /im and 
a dwelling time of thermal drawing « 100 sec, we can 
classify each materials combination by whether it falls in 
the r > Tdwciiing yellow region (radially stable) of Fig. [9] 
or in the t < Towelling white region (radially unstable). 

These materials combinations are presented in Fig. |10| 
The boundary line in red, which indicates the viscous 
combinations that satisfy t = Tdwciiing, divides the map 
into two areas. The shaded area above the boundary 
line is the region of potentially suitable materials combi- 
nations for fiber drawing (As2Se3-PES, AS2S3-PEI, Se- 
PSU) PUI [2T| : while the materials combinations below 
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the boundary line are unsuitable due to radial instability. 
Here, the only materials combination which seems to be 
definitely excluded in this regime due to radial instability 
is Se-PE. The polymer-polymer materials combination 
of PVDF-PC (polyvinylidene fluoride, PVDF, a piezo- 
electric polymer; polycarbonate, PC) is potentially suit- 
able for thermal drawing, and this possibility has been 
confirmed by recent experiments [39J. Moreover, since 
a high viscosity cladding improves stability (increase r), 
we predict that a wider variety of shell materials with 
low viscosity may possibly be employed in microstruc- 
tured fibers, such as the metals Sn and In |201l40j . These 
various available classes of metals, polymers and semicon- 
ductors expand the potential functionalities of devices in 
microstructured fibers. 



C. Stability of continuous filaments down to 
submicrometer /nanometer scale 

As shown in Fig. |2je-f), in some cases the initial 
shell breaks up (azimuthally) into long cylindrical fila- 
ments [211 [55]. These filaments themselves should be 
subject to the classic capillary instability and in prin- 
ciple should eventually break up into droplets. In our 
fiber-drawing experiments, however, no further instabil- 
ity is observed: the filaments are observed to remain con- 
tinuous and unbroken over at least cm length scales with 
diameters reaching submicrometer scales |211 122] . Since 
a cylindrical filament of one fiuid surrounded by another 
should be described exactly by the Tomotika model, it 
is important to know whether the observed stability is 
consistent with this model, or otherwise requires some 
additional physical mechanism (such as visco-elasticity) 
to explain. In this section, we find that the timescale 
of instability predicted by the Tomotika model exceeds 
the dwelling time of fiber drawing, making it unsurpris- 
ing that filament instability is not observed. In fact, we 
show that the instability time scale exceeds the dwelling 
time even under unrealistically conservative assumptions: 
that the filaments appear immediately in the fiber draw- 
ing, that the maximum instability growth rate is cumu- 
lative even though the lengthscale achieving maximum 
growth changes with the filament radius during drawing, 
and that the polymer viscosity is always at its minimum 
value (corresponding to the highest temperature point). 

An instability with time scale r corresponds to expo- 
nential growth of a fiuctuation amplitude e according to 
^ = e/r. If r is time varying, then the total ampli- 
tude growth is exp[J dt/T{t)]. Converting to dz — v{z)dt 
for a position-dependent axial fiow velocity v{z) during 
thermal drawing, we therefore obtain a total exponential 
growth factor: 



dz 



(9) 



/o v{z)t{z)' 

where z G [0, L] is axial position in the neck-down region 



with length (L = 6 cm). F ^ 1 corresponds to breakup, 
while F <C 1 corresponds to stability. 

In order to provide a conservative estimation of fila- 
ment stability, the capillary instability time is calculated 
from the fastest growth factor at each axial location (this 
is a very conservative estimate), and the polymer viscos- 
ity is set to be the minimum value (at the highest temper- 
ature) during thermal drawing. The capillary instability 
time scale is calculated based on the Tomotika model as 
follows |6], 



r(z) 



2^7^polymer(^!) ^) 



7{maXA (1 - .t2) $ [x, 77polymcr(z, i)/?7clad(2:, *)]} ' 

(10) 

The complex shape of neck-down profile is fitted from 
experiment can be approximately described by following 
formula, 



k = 



m 

R{L) 



-1, p = 2. (11) 



Due to the incompressibility of the viscous fiuid, the ve- 
locity of fiow scales inversely with area: 



v{0) ^ R^z)' 



(12) 



where v(0) = 4 x 10~''mm/sec is the preform velocity. 
Again by incompressibility, the filament radius (r) should 
scale as the fiber radius (R): 



r(0) 



R{oy 



(13) 



The temperature distribution during thermal drawing, fit 
from experiment, is found to be approximately parabolic. 



^ — ^max (-^max '^min) ^ 



2^-1 
L 



(14) 



In calculations, parameters for the typical As2Se3- 
PES fiber drawing are R(0) = 1 cm, s = 20, L = 6 cm, 
p = 2, T„,ax = 260°C, T„in-210°C, r(L) = 200 nm, 
'?poiymor — 10^ Pa-s. Fig. 11 (b)-(d) presents the corre- 



sponding position-dependent variables including radius, 
velocity, temperature and viscosity. Finally, we obtain 

F = 0.90. (15) 

This satisfies F < 1, but only barely — if this were an 
accurate estimate of the growth factor, instability might 
still be observed. However, the assumptions we made 
above were so conservative that the true growth factor 
must be much less than this, indicating the instability 
should not be observable during the dwelling time of fiber 
drawing. So, the observed filaments are consistent with 
the Tomotika model, although of course we cannot yet 
exclude the possibility that there are also additional ef- 
fects {e.g., elasticity) that further enhance stability. 
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Figure 11: Relevant parameters in the neck-down region dur- 
ing thermal drawing, (a) Photograph of neck-down region 
from preform to fiber, (b)-(e) for the calculated radius, ve- 
locity, temperature and viscosity. 



D. Favorability of azimuthal versus axial instability 



In experiments, we observed that thin films prefer- 
entially break up along the azimuthal direction rather 
than the axial direction. The discussion of the previ- 
ous Section [Vll C| suggests a simple geometrical explana- 
tion for such a preference, regardless of the details of the 
breakup mechanism. The key point is that any instability 
will have some characteristic wavelength A of maximum 
growth rate for small perturbations, and this A must be 
proportional to the characteristic feature size of the sys- 
tem, in this case the film thickness d. As the fiber is 
drawn, however, the thickness d and hence A decreases. 
Now, we consider what happens to an unstable pertur- 
bation that begins to grow at some wavelength Aq when 
the thickness is dg. If this is a perturbation along the ax- 
ial direction, then the fiber-draw process will stretch this 
perturbation to a longer wavelength, that will no longer 
correspond to the maximum-growth A (which is shrink- 
ing), and hence the growth will be damped. That is, the 
axial stretching competes with the layer shrinking, and 
will tend to suppress any axial breakup process. In con- 
trast, if Aq is an azimuthal perturbation, the draw-down 
process will shrink Aq along with the fiber cross-section 
at exactly the same rate that d and A shrink. There- 
fore, azimuthal instabilities are not suppressed by the 
draw process. This simple geometrical argument imme- 
diately predicts that the first observed instabilities will 
be azimuthal (although axial instabilities may still occur 
if the draw is sufficiently slow). 



VIII. CONCLUDING REMARKS 

In this paper, motivated by recent development in mi- 
crostructured optical fibers, we have explored capillary 
instability due to radial fiuctuations in a new geometry of 
concentric cylindrical shell by 2D numerical simulation, 
and applied its theoretical guidance to feature size and 
materials selections in the microstructured fibers during 
thermal drawing processes. Our results suggest several 
directions for future work. First, it would be desirable 
to extend the analytical theory of capillary instability 
in shells, which is currently available for equal viscos- 
ity only, to the more general case of unequal viscosi- 
ties — we have developed a very general theory in an 
appearing paper |25j . Second, we plan to extend our 
computation simulations to include 3D azimuthal fiuc- 
tuations together with radial fiuctuations; as argued in 
Section |VIID| we anticipate a general geometrical pref- 
erence for azimuthal breakup over axial breakup once 
the draw process is included. Third, there are many ad- 
ditional possible experiments that would be interesting 
to explore different aspects of these phenomena in more 
detail, such as employing different geometries {e.g., non- 
cylindrical), temperature-time profiles, or materials {e.g., 
Sn-PEI or Se-PE). Finally, by drawing more slowly so 
that axial breakup occurs, we expect that experiments 
should be able to obtain more diverse structures {e.g. ax- 
ial breakup into rings or complete breakup into droplets) 
that we hope to observe in the future. 
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Appendix 

1. Direct numerical simulation of concentric 
cylindrical shells 

In the simulation, the level set function is defined by 
a smoothed step function reinitialized at each time step 

{0 r<ri{z,t)~A, r > r2{z,t) + A 
1 ri{z,t) + A < r < r2{z,t) - A 
0.5 otherwise 

(16) 

where ri{z,t) and r2{z,t) are the radius of interfaces I 
and II of the coaxial cylinder, and A is the half thickness 
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of the discretized interface. The level set 4> = 1,0, and 0.5 
corresponds to the region of concentric cylindrical shell, 
of surrounding fluid, and interface, respectively. Contour 
n = {X\(p{X,t) — 0.5} tracks the interface I and II. 

A smoothed delta function is defined to project surface 
tension at interface. 



6 = 40(1 



1 |r-ri,2|<A 
otherwise 



A smoothed step function is introduced to create a 
smooth transition of the level-set function (f) from to 
1 across the interface, 



0.06 




0.00 



0.0 0.2 0.4 0.6 0.8 1.0 
Wavelength (27ira) 



1 + e{r-r^)/A 



1 + e{r^r,)/A 



(17) 



The boundary condition for the level set equation at 
the edge of the computational cell is 



Figure 12: Growth factor of instability as a function of per- 
turbation wavelength. Fast- and slow- modes occur at wave- 
lengths above their respective critical wavelengths A/ , As . In- 
set is a sketch of coa:xial cylinder with radius R = 2r and 
equal viscosities. 



ft ■ {~cV4> + (t>u) = 0, 



(18) 



where n and t are the normal and tangential vectors at 
the boundary. In addition, boundary conditions for the 
NS equations are 



t- ri\yu+ {\/uf]n =0. 



(19) 



In the simulation, time-stepping accuracy is controlled 
by an absolute and relative error tolerance (errabs and 
errrci) for each integration step. Let U be the solution 
vector at a given time step, E be the solver's estimated 
local error in U during this time step, and N be the 
number of degrees of freedom in the simulation. Then a 
time step is accepted if the following condition is satisfied. 



IS, 



\U^\ 



< 1. 



(20) 



A triangular finite-element mesh is generated, and 
second-order quadratic basis functions are used in the 
simulation. Parameters for Fig. [5] in the simulation are 
10^ kg/m^, 77 = 10^ Pa • s, 7 = 0.6 N/m, R = 120 
10-4, errabs = 10-^ D = Do = 10"" mVs. 



P 

/Ltm, errrei 



2. Linear theory of concentric cylindrical shells with 
equal viscosities 

A linear theory of capillary instability for a co-axial 
cylinder with equal viscosities is provided in the literature 
by Stone and Brenner [7J. The growth rate (a) for a wave 



vector k = 27r/A is a solution of the following quadratic 
equations. 



rrj 



"1 - [rkf] Kir, r) 



X < (7 - 



Rri 



[l-{Rkf] A{R, R) 



fc'^7i72 
rRrf 



'l-{rkf][l^{Rkf]k{r,Rf, (21) 



where r and R are the radii of the unperturbed inter- 
faces I and II, 7i and 72 are the interfacial tensions, and 
77 is viscosity. A(a, h), where a < &, is associated with the 
modified Bessel function, 



\f h\ f°° sJi{sa)Ji{sb) Id ,^r^ 

(22) 

For the case of 71 = 72 = 7, the growth rate has the 
following formula. 



a{X) 



1_ 
rjr 



*(A,i?/r), 



(23) 



where the growth factor of ^'(A, R/r) in Eq. 23 is a com- 
plicated function of instability wavelength [7]. The insta- 
bility time scale r ^ ~ rjr/^ is scaled with radius. 
For the case of i? = 2r, this growth factor is calculated 
in Fig. [12] A positive growth factor indicates a pos- 
itive growth rate (cr > 0), for which any perturbation 
is exponentially amplified with time. Instability occurs 
at long wavelengths above a certain critical wavelength. 
Two critical wavelengths exist for the co-axial cylinder 
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Figure 13: Temperature-dependent viscosity for various 
chalcogenide glasses. Typical temperature during fiber draw- 
ing for glass Se, As2Se3, AS2S3 is around 220, 260, 300 °C with 
the corresponding viscosities of 10,10^,10^ Pa • s, respec- 
tively. 

shell. One is a short critical wavelength A/ = 27rr for a 
faster-growth mode (red line). The other is a long crit- 
ical wavelength = 27ri? for slower-growth mode (blue 
line). In the numerical simulation, the wavelength is cho- 
sen between these two wavelengths (A/ < A < As) , and 
fast modes dominate. From the simulation parameters 
7y sa 10^ Pa • s and 7/r ~ 10^ Pa • s, together with the 
wavelength 27rr/A ~ 0.47 corresponding to a growth fac- 
tor \E'(A) « 0.03, the linear theory predicts an instability 
time scale r = cr"^ = {^^{X, R/r))'^ w 334 sec. 



3. Viscous materials during tliermal drawing 

Our chosen materials include chalcogenide glasses (Se, 
As2Se3, and AS2S3) and thermoplastic polymers (PES, 
PEI, and PSU). The viscosity of chalcogenide glass- 
forming melts depends on temperature and is calculated 
from an empirical Arrhenius formula|42| . 



logr; = logr;o + ~ 1, (24) 

where R is the ideal gas constant, T is the tempera- 
ture in Kelvin, and r] is viscosity in Pa • s. The param- 
eters of logyyo, C, and D for our materials are listed be- 
low: -2.0,6651,770.82 for Se, -3.09,18877.8,875.56 for 
As2Se3, and -3.62, 33744, 650.8 for AsaSaPSl. These vis- 
cosities over a wide temperature range are plotted in Fig. 
[TS] The typical temperature during a fiber drawing for 
Se,As2Se3, or AS2S3 films is around 220,260, or 300 °C, 
respectively, with the corresponding viscosities of 10, 10^, 
or 10^ Pa • s, respectively. 

4. Diffusion term in level-set equation 

To ensure numerical stability in the simulation, an ar- 
tificial diffusion term proportional to a small parameter 
D is added to level-set Eq. [5] as follows HI] , 



0t-f u. V(/) = i:>v-{[i-(/)(i-(?;.)]V(/)}. (25) 
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